package com.example.analysis.utils;

import com.example.analysis.entity.Quest;
import com.example.analysis.exception.BusinessException;
import com.example.analysis.vo.ValidityFactorLoading;
import com.example.analysis.vo.ValidityVO;
import com.example.analysis.vo.VarianceExplained;
import org.apache.commons.math3.distribution.ChiSquaredDistribution;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Comparator;
import java.util.List;

public class FactorAnalysis {
    public static ValidityVO analysis(double[][] data, List<Quest> questList, int factorNum) throws BusinessException {
        double[][] normalizedData = ComCalcUtils.zScore(data);
        // 1. 检验
        // 1.1 计算KMO值
        double[][] correlationMatrix = ComCalcUtils.calculateCorrelationMatrix(normalizedData);
        double[][] partialCorrelationMatrix = ComCalcUtils.calculatePartialCorrelationMatrix(correlationMatrix);
        double kmo = calculateKMO(correlationMatrix, partialCorrelationMatrix);
        System.out.println("kmo = " + kmo);
        // 1.2 Bartlett 球形度检验
        double[] bartlettTest = bartlettTest(normalizedData);
        if (bartlettTest[1] >= 0.05) {
            System.out.println("不适合因子分析");
            return new ValidityVO(kmo, bartlettTest[0], bartlettTest[1]);
        }
        // 2. 计算协方差
        double[][] covarianceMatrix = ComCalcUtils.calculateCovarianceMatrix(normalizedData);
        // 3. 求特征值和特征向量
        EigenDecomposition eigenDecomposition = calcVectorsAndValues(covarianceMatrix);
        //   3.1 特征值
        double[] eigenvalues = eigenDecomposition.getRealEigenvalues();
        RealMatrix eigenVectors = eigenDecomposition.getV();
        //   3.2 按特征值降序排序
        Integer[] indices = sortEigenvaluesDescending(eigenvalues);
        eigenvalues = reorderArray(eigenvalues, indices);
        eigenVectors = reorderMatrixColumns(eigenVectors, indices);
        // 4. 计算因子数量
        if (factorNum == 0)
            factorNum = calcFactorNum(eigenvalues);
        // 5. 计算指标
        //   5.1 计算因子荷载矩阵
        double[][] factorLoadings = extractFactorLoadings(eigenVectors, eigenvalues, covarianceMatrix.length, factorNum);
        //   5.2 计算方差解释率和累计方差解释率
        double[][] explained = calculateVarianceExplained(eigenvalues);
        // 6. 因子旋转
        double[][] rotatedLoadings = rotateFactors(factorLoadings);
        // 7. 计算旋转后的方差解释率和累计方差解释率
        double[][] rotatedVarianceExplained = calculateRotatedVarianceExplained(rotatedLoadings);

        // 封装
        // 方差解释率
        List<VarianceExplained> varianceExplainedList = new ArrayList<>();
        for (int i = 0; i < eigenvalues.length; i++) {
            varianceExplainedList.add(new VarianceExplained(eigenvalues[i], explained[0][i], explained[1][i]));
        }
        // 旋转后方差解释率
        List<VarianceExplained> rotatedVarianceExplainedList = new ArrayList<>();
        for (int i = 0; i < rotatedVarianceExplained[0].length; i++) {
            rotatedVarianceExplainedList.add(new VarianceExplained(rotatedVarianceExplained[0][i], rotatedVarianceExplained[1][i], rotatedVarianceExplained[2][i]));
        }
        // 因子荷载
        List<ValidityFactorLoading> factorLoadingList = new ArrayList<>();
        for (int i = 0; i < rotatedLoadings.length; i++) {
            factorLoadingList.add(new ValidityFactorLoading(questList.get(i).getQuestEvaluate(), rotatedLoadings[i]));
        }
        return new ValidityVO(kmo, bartlettTest[0], bartlettTest[1], varianceExplainedList, rotatedVarianceExplainedList, factorLoadingList);
    }

    // Bartlett 球形度检验
    private static double[] bartlettTest(double[][] data) throws BusinessException {
        // 计算相关系数矩阵
        double[][] correlationMatrix = ComCalcUtils.calculateCorrelationMatrix(data);
        // 转成矩阵
        Array2DRowRealMatrix corrMatrix = new Array2DRowRealMatrix(correlationMatrix);
        // 计算相关矩阵的行列式（需处理奇异矩阵）
        LUDecomposition luDecomp = new LUDecomposition(corrMatrix);
        double det = luDecomp.getDeterminant();
        if (det <= 0) {
             throw new BusinessException("相关矩阵不可逆，无法计算Bartlett检验");
        }
        // Bartlett统计量公式
        int n = data[0].length;      // 样本量
        int p = data.length;   // 变量数
        double m = n - 1 - (2 * p + 5) / 6.0;
        double chi2 = -m * Math.log(det);
        int df = p * (p - 1) / 2;  // 自由度
        // 计算p值
        ChiSquaredDistribution chiDist = new ChiSquaredDistribution(df);
        double pValue = 1 - chiDist.cumulativeProbability(chi2);
        System.out.println("Bartlett 球形度检验：统计量" + chi2 + ", pValue = " + pValue);
        return new double[]{chi2, pValue};
    }

    // 计算 KMO 值
    public static double calculateKMO(double[][] correlationMatrix, double[][] partialCorrelationMatrix) {
        int p = correlationMatrix.length; // 题项数量
        double sumRijSquared = 0;
        double sumPijSquared = 0;

        for (int i = 0; i < p; i++) {
            for (int j = 0; j < p; j++) {
                if (i != j) {
                    sumRijSquared += correlationMatrix[i][j] * correlationMatrix[i][j];
                    sumPijSquared += partialCorrelationMatrix[i][j] * partialCorrelationMatrix[i][j];
                }
            }
        }
        return sumRijSquared / (sumRijSquared + sumPijSquared);
    }


    // ---求特征值和特征向量
    private static EigenDecomposition calcVectorsAndValues(double[][] covarianceMatrix) {
        RealMatrix matrix = new Array2DRowRealMatrix(covarianceMatrix);
        return new EigenDecomposition(matrix);
    }

    // 计算方差解释率和累积方差解释率
    public static double[][] calculateVarianceExplained(double[] eigenValues) {
        double totalVariance = 0;
        for (double eigenValue : eigenValues) {
            totalVariance += eigenValue;
        }

        int numFactors = eigenValues.length;
        double[] varianceExplained = new double[numFactors];
        double[] cumulativeVarianceExplained = new double[numFactors];

        System.out.println("因子\t特征值\t方差解释率\t累积方差解释率");
        for (int i = 0; i < numFactors; i++) {
            varianceExplained[i] = eigenValues[i] / totalVariance;
            if (i == 0) {
                cumulativeVarianceExplained[i] = varianceExplained[i];
            } else {
                cumulativeVarianceExplained[i] = cumulativeVarianceExplained[i - 1] + varianceExplained[i];
            }
            System.out.printf("%d\t%.4f\t%.4f\t%.4f\n", i + 1, eigenValues[i], varianceExplained[i], cumulativeVarianceExplained[i]);
        }
        return new double[][]{varianceExplained, cumulativeVarianceExplained};
    }

    // 计算旋转后的方差解释率和累积方差解释率
    public static double[][] calculateRotatedVarianceExplained(double[][] rotatedLoadings) {
        int p = rotatedLoadings.length;// 题项数量，也是旋转后总方差（标准化后）
        int m = rotatedLoadings[0].length;// 因子数量
        double[] factorVariances = new double[m];

        for (int j = 0; j < m; j++) {
            double sum = 0;
            for (double[] rotatedLoading : rotatedLoadings) {
                sum += rotatedLoading[j] * rotatedLoading[j];
            }
            factorVariances[j] = sum;
        }

        double[] varianceExplained = new double[m];
        double[] cumulativeVarianceExplained = new double[m];

        for (int i = 0; i < m; i++) {
            varianceExplained[i] = factorVariances[i] / p;
            if (i == 0) {
                cumulativeVarianceExplained[i] = varianceExplained[i];
            } else {
                cumulativeVarianceExplained[i] = cumulativeVarianceExplained[i - 1] + varianceExplained[i];
            }
        }

        System.out.println("旋转后因子\t特征值\t方差解释率\t累积方差解释率");
        for (int i = 0; i < m; i++) {
            System.out.printf("%d\t%.4f\t%.4f\t%.4f\n", i + 1, factorVariances[i], varianceExplained[i], cumulativeVarianceExplained[i]);
        }
        return new double[][]{factorVariances, varianceExplained, cumulativeVarianceExplained};
    }

    // 确定因子数量。特征值 > 1的数量
    public static int calcFactorNum(double[] eigenValues) {
        int factorNum = 0;
        for (double v : eigenValues) {
            if (v > 1)
                factorNum++;
        }
        return factorNum;
    }

    // ----方差最大旋转
    public static double[][] rotateFactors(double[][] factorLoadings) {
        int p = factorLoadings.length; // 题项数量
        int m = factorLoadings[0].length; // 因子数量
        double[][] rotatedLoadings = new double[p][m];
        for (int i = 0; i < p; i++) {
            System.arraycopy(factorLoadings[i], 0, rotatedLoadings[i], 0, m);
        }
        int maxIterations = 100;
        double tolerance = 1e-6;
        double delta = Double.MAX_VALUE;
        int iteration = 0;
        while (delta > tolerance && iteration < maxIterations) {
            double oldDelta = delta;
            delta = 0;
            for (int k = 0; k < m - 1; k++) {
                for (int l = k + 1; l < m; l++) {
                    double sum1 = 0, sum2 = 0, sum3 = 0, sum4 = 0;
                    for (int i = 0; i < p; i++) {
                        double hk = rotatedLoadings[i][k];
                        double hl = rotatedLoadings[i][l];
                        double hk2 = hk * hk;
                        double hl2 = hl * hl;
                        double hk4 = hk2 * hk2;
                        double hl4 = hl2 * hl2;
                        sum1 += (hk2 - hl2) * hk * hl;
                        sum2 += (hk4 - hl4);
                        sum3 += (3 * hk2 - hl2) * hk2;
                        sum4 += (3 * hl2 - hk2) * hl2;
                    }
                    // 计算旋转角度
                    double gamma = (sum3 - sum4) / (4 * sum1);
                    double sinPhi = Math.sin(2 * gamma);
                    double cosPhi = Math.cos(2 * gamma);
                    // 更新。旋转因子对 k,l
                    for (int i = 0; i < p; i++) {
                        double hk = rotatedLoadings[i][k];
                        double hl = rotatedLoadings[i][l];
                        rotatedLoadings[i][k] = hk * cosPhi + hl * sinPhi;
                        rotatedLoadings[i][l] = -hk * sinPhi + hl * cosPhi;
                    }
                    delta += Math.abs(gamma);
                }
            }
            if (oldDelta == delta) {
                break;
            }
            iteration++;
        }
        System.out.println("旋转后因子荷载");
        for (int i = 0; i < p; i++) {
            System.out.printf("因子%d\t", (i + 1));
            for (int j = 0; j < m; j++) {
                System.out.print(rotatedLoadings[i][j] + "\t");
            }
            System.out.println();
        }
        return rotatedLoadings;
    }

    // ----提取因子载荷矩阵
    // p: 题项，numFactors: 要提取的因子数量
    // 按特征值降序排序并返回索引
    private static Integer[] sortEigenvaluesDescending(double[] eigenvalues) {
        Integer[] indices = new Integer[eigenvalues.length];
        for (int i = 0; i < indices.length; i++) indices[i] = i;
        Arrays.sort(indices, Comparator.comparingDouble((Integer i) -> -eigenvalues[i]));
        return indices;
    }

    // 根据索引重新排序数组
    private static double[] reorderArray(double[] arr, Integer[] indices) {
        double[] reordered = new double[arr.length];
        for (int i = 0; i < arr.length; i++) {
            reordered[i] = arr[indices[i]];
        }
        return reordered;
    }

    // 根据索引重新排序矩阵列
    private static RealMatrix reorderMatrixColumns(RealMatrix matrix, Integer[] indices) {
        RealMatrix reordered = new Array2DRowRealMatrix(matrix.getRowDimension(), matrix.getColumnDimension());
        for (int i = 0; i < indices.length; i++) {
            reordered.setColumn(i, matrix.getColumn(indices[i]));
        }
        return reordered;
    }

    // 荷载矩阵
    public static double[][] extractFactorLoadings(RealMatrix eigenVectors, double[] eigenValues, int p, int numFactors) {
        double[][] factorLoadings = new double[p][numFactors];
        System.out.println("因子载荷");
        for (int i = 0; i < p; i++) {
            for (int j = 0; j < numFactors; j++) {
                factorLoadings[i][j] = eigenVectors.getEntry(i, j) * Math.sqrt(eigenValues[j]);
                System.out.print(factorLoadings[i][j] + ",");
            }
            System.out.println();
        }

        return factorLoadings;
    }
}
